import math, copy, numpy as np, scipy.stats as ss, warnings warnings.filterwarnings("ignore") def power_pooled(n,C,alpha,gamma,tau,sigma,pi,f): # Critical t-values t_alpha = ss.t.ppf(1 - alpha/2,n*C-3) # critical value for Type I error t_gamma = ss.t.ppf(gamma,n*C-3) # critical value for Type II error # Initialize distribution statistics mu_ind = np.zeros(len(pi)) eta_ind = np.zeros(len(pi)) # Calculate distribution statistics for i in range(len(pi)): mu_ind[i] = pi[i] * f[i] eta_ind[i] = pi[i]**2 * f[i] mu = sum(mu_ind) eta = sum(eta_ind) # eta=E[pi^2] here; eta=eta-mu^2 in the paper. psi = 0 if(pi[0] == 0): psi = f[0] muS = 1 - mu - psi etaT = (eta-mu**2)/(1-psi)-(psi/(1-psi)**2)*mu**2 # Variance multipliers varN = tau + sigma varCo = (n-1) * tau ########### Minimum Detectable Effects ######### # In the case of a non-trivial design Var = 1/(n*C)*(varCo*(1/(psi*(1-psi)) + (1-psi)/(mu **2)*etaT) + varN*(psi+mu )/(mu *psi)) VarS = 1/(n*C)*(varCo*(1/(psi*(1-psi)) + (1-psi)/(muS**2)*etaT) + varN*(psi+muS)/(muS*psi)) MDE_T = (t_alpha + t_gamma)*Var**0.5 MDE_S = (t_alpha + t_gamma)*VarS**0.5 # In the case of a pure control Var_T = 1/(n*C)*(varCo*(eta-mu**2)/(mu**2*(1-mu)**2) + varN/(mu*(1-mu))) MDE_Tonly = (t_alpha+t_gamma)*Var_T**0.5 # In the case of a pure treatment Var_S = 1/(n*C)*(varCo*(eta-(1-mu)**2)/(mu**2*(1-mu)**2) + varN/(mu*(1-mu))) MDE_Sonly = (t_alpha+t_gamma)*Var_S**0.5 return {"MDE_T":MDE_T, "MDE_S":MDE_S, "MDE_Tonly":MDE_Tonly} def power_slope(n,C,alpha,gamma,tau,sigma,pi,f,j,k): # Critical t-values t_alpha = ss.t.ppf(1 - alpha/2,n*C-3) # critical value for Type I error t_gamma = ss.t.ppf(gamma,n*C-3) # critical value for Type II error # Initialize distribution statistics mu_ind = np.zeros(len(pi)) p_ind = np.zeros(len(pi)) # Calculate distribution statistics for i in range(len(pi)): mu_ind[i] = pi[i] * f[i] p_ind[i] = (1-pi[i]) * f[i] # Variance multipliers varN = tau + sigma varCo = (n-1) * tau ######### Minimum Detectable Effects ######### Var_T = (varCo*(1/f[j]+1/f[k])+varN*(1/mu_ind[j]+1/mu_ind[k]))/(n*C) Var_S = (varCo*(1/f[j]+1/f[k])+varN*(1/ p_ind[j]+1/ p_ind[k]))/(n*C) MDSE_T = ((t_alpha+t_gamma)/(pi[k]-pi[j]))*Var_T**0.5 MDSE_S = ((t_alpha+t_gamma)/(pi[k]-pi[j]))*Var_S**0.5 return {"MDSE_T":MDSE_T, "MDSE_S":MDSE_S} def power_ind(n,C,alpha,gamma,tau,sigma,pi,f,p): # Critical t-values t_alpha = ss.t.ppf(1 - alpha/2,n*C-3) # critical value for Type I error t_gamma = ss.t.ppf(gamma,n*C-3) # critical value for Type II error # Initialize distribution statistics mu_ind = np.zeros(len(pi)) p_ind = np/zeros(len(pi)) # Calculate distribution statistics for i in range(len(pi)): mu_ind[i] = pi[i] * f[i] p_ind[i] = (1-pi[i]) * f[i] psi = 0 if(pi[0] == 0): psi = f[0] # Variance multipliers varN = tau + sigma varCo = (n-1)*tau varC = n*tau + sigma ########### Minimum Detectable Effects ######### MDE_ind_T = (t_alpha+t_gamma)*(1/(n*C)*(varCo*(1/f[p]+1/psi)+varN*(1/mu_ind[p]+1/psi)))**0.5 MDE_ind_S = (t_alpha+t_gamma)*(1/(n*C)*(varCo*(1/f[p]+1/psi)+varN*(1/p_ind[p] +1/psi)))**0.5 return {"MDE_ind_T":MDE_ind_T, "MDE_ind_S":MDE_ind_S} def Compass_Search(f, admissable, x0, delta = 1.0, tol=1e-5): # INPUTS: # f is a function on domain(x) with list output containing at least $obj to minimize. Ancillary output is optional. # admissable is a function on domain(x), and returns a boolean if that value is allowed. # x0 is the original guess for the optimal x. # delta is the initial sizes of changes made in guesses for x. # tol is the maximum error tolerated in the final answer for. # OUTPUTS: # x is the minimizer of f. # y is the minimum $obj value of f(x). # steps is the number of steps is took to complete the optimization, up to a multiplicative constant. def f_obj(x): return f(x)["obj"] ##xs = [] # optional object to store working values of x. d = len(x0) D = np.vstack([np.diag([delta]*d), np.diag([-delta]*d)]) steps = 0 x = x0 while D[0,0] > tol: steps += 1 admissable_rows = np.where(np.apply_along_axis(admissable, 1, x + D))[0] # check which rows are admissable. if len(admissable_rows) == 0: D /= 2 else: explore_values = np.apply_along_axis(f_obj, 1, x + D[admissable_rows,:]) min_index = int(np.where(explore_values == min(explore_values))[0][0]) no_improvement = explore_values[min_index] >= f(x)["obj"] # if all values are not an improvement, if no_improvement: D = D/2 else: x = x + D[admissable_rows[min_index],] ##xs.append(x) return {"x":x, "y":f(x), "steps":steps*d}#, "xs":xs def index_min(x): # index of minimum of matrix x. This can be used for brute search, which is commented out here. # index of minimum of matrix x i = 0 # row j = 0 # column rows = x.shape[0] columns = x.shape[1] for r in range(rows): for c in range(columns): if x[r,c] < x[i,j]: i = r j = c return i, j # ################################################# # # # Replication Code for Baird et al. 2016 # # # # ################################################# # ################################################# # Table 1, Column 1 # # ################################################# n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0 # variance of cluster error (tau^2) sigma = 1 # variance of individual error (sigma^2) pi = [0,1] # cluster saturations f = [1/2,1/2] # fraction of clusters in each saturation, respectively X = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) print(X["MDE_Tonly"]) # ################################################# # Table 1, Column 2 # # ################################################# n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) pi = [0,1] # cluster saturations f = [1/2,1/2] # fraction of clusters in each saturation, respectively X = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) print(X["MDE_Tonly"]) # ################################################# # Table 1, Column 3 # # ################################################# n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 1.0 # variance of cluster error (tau^2) sigma = 0.0 # variance of individual error (sigma^2) pi = [0,1] # cluster saturations f = [1/2,1/2] # fraction of clusters in each saturation, respectively X = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) print(X["MDE_Tonly"]) # ################################################# # Table 1, Column 4 # # ################################################# # Partial Population design # Half of clusters pure control, half partially treated, ICC = 0.0 n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0 # variance of cluster error (tau^2) sigma = 1 # variance of individual error (sigma^2) pi = [0,0.5] # cluster saturations theta = 1/2 # proportion of optimization weight on MDE_T # BRUTE FORCE SEARCH ##psi = np.arange(0.01,.99,.001) ##X = np.zeros((len(psi), 3)) ##for j in range(len(psi)): ## f = [psi[j],1-psi[j]] ## pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## X[j,] = [pooled[i] for i in ["MDE_T","MDE_S","MDE_Tonly"]] ## ##obj = X[:,:2] * [theta, 1-theta] # objective function to maximize ##x = index_min(obj) # index of psi that minimizes objective ##print(psi[x[0]]) # optimal proportion of pure control clusters ##print(1 - psi[x[0]]) # optimal proportion of partial treatment clusters ##print(X[x[0],0]) # MDE_T at optimal control size ##print(X[x[0],1]) # MDE_S at optimal control size # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0]<1 def objective(x): # x contains the first saturation level. f = [x[0],1-x[0]] value = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) objective_function = theta*value["MDE_T"] + (1-theta)*value["MDE_S"] # objective function to maximize: a linear combination of MDE_T and MDE_S. return {"obj":objective_function, "anc":{"optMDE":value}} opt = Compass_Search(objective, admissable, x0 = np.array([.5])) print(opt["x"][0]) # optimal proportion of pure control clusters print(1-opt["x"][0]) # optimal proportion of partial treatment clusters print(opt["y"]["anc"]["optMDE"]["MDE_T"]) # MDE_T at optimal control size print(opt["y"]["anc"]["optMDE"]["MDE_S"]) # MDE_S at optimal control size # ################################################# # Table 1, Column 5 # # ################################################# # Partial Population design # Half of clusters pure control, half partially treated, ICC = 0.1 n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) pi = [0,0.5] # cluster saturations theta = 1/2 # proportion of optimization weight on MDE_T # BRUTE FORCE SEARCH ##psi = np.arange(0.01,.99,.001) ##X = np.zeros((len(psi), 3)) ##for j in range(len(psi)): ## f = [psi[j],1-psi[j]] ## pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## X[j,] = [pooled[i] for i in ["MDE_T","MDE_S","MDE_Tonly"]] ##obj = X[:,:2] * [theta, 1-theta] # objective function to maximize ##x = index_min(obj) # index of psi that minimizes objective ##print(psi[x[0]]) # optimal proportion of pure control clusters ##print(1-psi[x[0]]) # optimal proportion of partial treatment clusters ##print(X[x[0],0]) # MDE_T at optimal control size ##print(X[x[0],1]) # MDE_S at optimal control size # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0]<1 def objective(x): # x contains the first saturation level. f = [x[0],1-x[0]] value = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) objective_function = theta*value["MDE_T"] + (1-theta)*value["MDE_S"] # objective function to maximize: a linear combination of MDE_T and MDE_S. return {"obj":objective_function, "anc":{"optMDE":value}} opt = Compass_Search(objective, admissable, x0 = np.array([.5])) print(opt["x"][0]) # optimal proportion of pure control clusters print(1-opt["x"][0]) # optimal proportion of partial treatment clusters print(opt["y"]["anc"]["optMDE"]["MDE_T"]) # MDE_T at optimal control size print(opt["y"]["anc"]["optMDE"]["MDE_S"]) # MDE_S at optimal control size # ################################################# # Table 1, Column 6 # # ################################################# # Partial Population design # Half of clusters pure control, half partially treated, ICC = 1.0 n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 1.0 # variance of cluster error (tau^2) sigma = 0.0 # variance of individual error (sigma^2) pi = [0,0.5] # cluster saturations theta = 1/2 # proportion of optimization weight on MDE_T # BRUTE FORCE SEARCH ##psi = np.arange(0.01,.99,.001) ##X = np.zeros((len(psi), 3)) ##for j in range(len(psi)): ## f = [psi[j],1-psi[j]] ## pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## X[j,] = [pooled[i] for i in ["MDE_T","MDE_S","MDE_Tonly"]] ##obj = X[:,:2] * [theta, 1-theta] # objective function to maximize ##x = index_min(obj) # index of psi that minimizes objective ##print(psi[x[0]]) # optimal proportion of pure control clusters ##print(1-psi[x[0]]) # optimal proportion of partial treatment clusters ##print(X[x[0],0]) # MDE_T at optimal control size ##print(X[x[0],1]) # MDE_S at optimal control size # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0]<1 def objective(x): # x contains the first saturation level. f = [x[0],1-x[0]] value = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) objective_function = theta*value["MDE_T"] + (1-theta)*value["MDE_S"] # objective function to maximize: a linear combination of MDE_T and MDE_S. return {"obj":objective_function, "anc":{"optMDE":value}} opt = Compass_Search(objective, admissable, x0 = np.array([.5])) print(opt["x"][0]) # optimal proportion of pure control clusters print(1-opt["x"][0]) # optimal proportion of partial treatment clusters print(opt["y"]["anc"]["optMDE"]["MDE_T"]) # MDE_T at optimal control size print(opt["y"]["anc"]["optMDE"]["MDE_S"]) # MDE_S at optimal control size # ################################################# # Table 1, Column 7-8 # # ################################################# #### Table 1, Column 7 # Partial Population design # Some clusters pure control, the rest partially treated, ICC = 0.1 # Subject to twice as much weight on MDE_S compared to MDE_T n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) theta = 1/3 # proportion of optimization weight on MDE_T. # BRUTE FORCE SEARCH ##x = np.arange(0.01,.99,.005) # values for optimization ##MDE_T = np.zeros((len(x),len(x))) ##MDE_S = copy.deepcopy(MDE_T) ##for i in range(len(x)): ## pi = [0,x[i]] ## for j in range(len(x)): ## f = [x[j],1-x[j]] ## pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## X = [pooled[k] for k in ["MDE_T","MDE_S","MDE_Tonly"]] ## MDE_T[i,j] = X[0] # MDE_T ## MDE_S[i,j] = X[1] # MDE_S ## ## ##obj = theta*MDE_T + (1-theta)*MDE_S # objective function to maximize ##ind = index_min(obj) # index of psi that minimizes MDE_T ##print(x[ind[0]]) # optimal saturation of treatment ##print(x[ind[1]]) # optimal proportion of pure control clusters ##print(1-x[ind[1]]) # optimal proportion of partial treatment clusters ##print(MDE_T[ind[0],ind[1]]) # MDE_T at optimal control size ##print(MDE_S[ind[0],ind[1]]) # MDE_T at optimal control size # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0] < 1 and 0 < x[1] and x[1] < 1 def objective(x):# x contains the second saturation level, and the number of people in the pure control group. pi = [0, x[0]] f = [x[1], 1-x[1]] value = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) objective_function = theta*value["MDE_T"] + (1-theta)*value["MDE_S"] # objective function to maximize: a linear combination of MDE_T and MDE_S. return {"obj":objective_function, "anc":{"optMDE":value}} opt = Compass_Search(objective, admissable, x0 = np.array([.5,.5])) print(opt["x"][0]) # optimal saturation of treatment print(opt["x"][1]) # optimal proportion of pure control clusters print(1-opt["x"][1]) # optimal proportion of partial treatment clusters print(opt["y"]["anc"]["optMDE"]["MDE_T"]) # MDE_T at optimal control size print(opt["y"]["anc"]["optMDE"]["MDE_S"]) # MDE_T at optimal control size #### Table 1, Column 8 # Partial Population design # Some clusters pure control, the rest partially treated, ICC = 0.1 # Subject to MDE_T <= 0.25 # BRUTE FORCE SEARCH ##MDE_S_aug = copy.deepcopy(MDE_S) ##for i in range(len(x)): ## for j in range(len(x)): ## if(MDE_T[i,j] > 0.25): ## MDE_S_aug[i,j] = np.inf ## ##obj = MDE_S_aug # objective function to maximize ##ind = index_min(obj) # index of psi that minimizes MDE_S ##print(x[ind[0]]) # optimal saturation of treatment ##print(x[ind[1]]) # optimal proportion of pure control clusters ##print(1 - x[ind[1]]) # optimal proportion of partial treatment clusters ##print(MDE_T[ind[0], ind[1]]) # MDE_T at optimal control size ##print(MDE_S[ind[0], ind[1]]) # MDE_T at optimal control size # COMPASS SEARCH def minimum_admissable(x): return 0 < x and x < 1 def minimum_objective(theta): def admissable(x): return 0 < x[0] and x[0] < 1 and 0 < x[1] and x[1]<1 def objective(x):# x contains the second saturation level, and the number of people in the pure control group. pi = [0, x[0]] f = [x[1], 1-x[1]] pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) objective_function = theta*pooled["MDE_T"] + (1-theta)*pooled["MDE_S"] return {"obj":objective_function, "anc":{"x":x, "pooled":pooled}} opt = Compass_Search(objective, admissable, x0 = np.array([0.5, 0.5])) return {"obj": abs(opt["y"]["anc"]["pooled"]["MDE_T"]-0.25), "anc": {"opt":opt}} theta = Compass_Search(minimum_objective, minimum_admissable, x0=np.array([.5])) print(theta["y"]["anc"]["opt"]["x"][0]) # optimal saturation of treatment print(theta["y"]["anc"]["opt"]["x"][1]) # optimal proportion of pure control clusters print(1-theta["y"]["anc"]["opt"]["x"][1]) # optimal proportion of partial treatment clusters print(theta["y"]["anc"]["opt"]["y"]["anc"]["pooled"]["MDE_T"]) # MDE_T at optimal control size print(theta["y"]["anc"]["opt"]["y"]["anc"]["pooled"]["MDE_S"]) # MDE_S at optimal control size # ################################################# # Table 2, Column 1 # # ################################################# # Random saturation design (two partially treated clusters), ICC = 0.0 # Optimal cluster saturations and allocation n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.0 # variance of cluster error (tau^2) sigma = 1.0 # variance of individual error (sigma^2) theta = 1/2 # proportion of optimization weight on MDSE_T # BRUTE FORCE SEARCH ##x = np.arange(0.01,.5+.005,.01) ##MDSE_T = np.zeros((len(x),len(x))) ##MDSE_S = copy.deepcopy(MDSE_T) ##for i in range(len(x)): ## pi = [x[i],1-x[i]] ## for j in range(len(x)): ## f = [x[j],1-x[j]] ## slope = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) ## X = [slope[k] for k in ["MDSE_T","MDSE_S"]] ## MDSE_T[i,j] = copy.deepcopy(X[0]) ## MDSE_S[i,j] = copy.deepcopy(X[1]) ## ##obj = theta*MDSE_T + (1-theta)*MDSE_S; # objective function to maximize ##ind = index_min(obj); # index of psi that minimizes objective ##print(x[ind[0]]) # optimal saturation for lower saturation ##print(1 - x[ind[0]]) # optimal saturation for higher saturation ##print(x[ind[1]]) # optimal fraction of clusters in lower saturation ##print(1 - x[ind[1]]) # optimal fraction of clusters in higher saturation ##print(MDSE_T[ind[0], ind[1]]) # MDSE_T at optimum ##print(MDSE_S[ind[0], ind[1]]) # MDSE_S at optimum # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0] < .5 and 0 < x[1] and x[1] < .5 def objective(x): # x contains the first saturation level, and the number of people in the least saturated group. pi = [x[0],1-x[0]] f = [x[1],1-x[1]] value = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) objective_function = theta*value["MDSE_T"] + theta*value["MDSE_S"] return {"obj":objective_function, "anc":{"opt_MDSE":value}} opt = Compass_Search(objective, admissable, x0=np.array([.25,.25])) print(opt["x"][0]) # optimal saturation for lower saturation print(1-opt["x"][0]) # optimal saturation for higher saturation print(opt["x"][1]) # optimal fraction of clusters in lower saturation print(1-opt["x"][1]) # optimal fraction of clusters in higher saturation print(opt["y"]["anc"]["opt_MDSE"]["MDSE_T"]) # MDSE_T at optimum print(opt["y"]["anc"]["opt_MDSE"]["MDSE_S"]) # MDSE_S at optimum # ################################################# # Table 2, Column 2 # # ################################################# # Random saturation design (two partially treated clusters), ICC = 0.0 # Optimal cluster saturations and allocation n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) # BRUTE FORCE SEARCH ##x = np.arange(0.01,.5+.005,.01) ##MDSE_T = np.zeros((len(x),len(x))) ##MDSE_S = copy.deepcopy(MDSE_T) ##for i in range(len(x)): ## pi = [x[i],1-x[i]] ## for j in range(len(x)): ## f = [x[j],1-x[j]] ## slope = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) ## X = [slope[k] for k in ["MDSE_T","MDSE_S"]] ## MDSE_T[i,j] = copy.deepcopy(X[0]) ## MDSE_S[i,j] = copy.deepcopy(X[1]) ## ## ##theta = 1/2 # fraction of weight put onto MDSE_T ##obj = theta*MDSE_T + (1-theta)*MDSE_S; # objective function to maximize ##ind = index_min(obj); # index of psi that minimizes objective ##print(x[ind[0]]) # optimal saturation for lower saturation ##print(1 - x[ind[0]]) # optimal saturation for higher saturation ##print(x[ind[1]]) # optimal fraction of clusters in lower saturation ##print(1 - x[ind[1]]) # optimal fraction of clusters in higher saturation ##print(MDSE_T[ind[0], ind[1]]) # MDSE_T at optimum ##print(MDSE_S[ind[0], ind[1]]) # MDSE_S at optimum # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0] < .5 and 0 < x[1] and x[1] < .5 def objective(x): # x contains the first saturation level, and the number of people in the least saturated group. pi = [x[0],1-x[0]] f = [x[1],1-x[1]] value = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) objective_function = theta*value["MDSE_T"] + theta*value["MDSE_S"] return {"obj":objective_function, "anc":{"opt_MDSE":value}} opt = Compass_Search(objective, admissable, x0=np.array([.25,.25])) print(opt["x"][0]) # optimal saturation for lower saturation print(1-opt["x"][0]) # optimal saturation for higher saturation print(opt["x"][1]) # optimal fraction of clusters in lower saturation print(1-opt["x"][1]) # optimal fraction of clusters in higher saturation print(opt["y"]["anc"]["opt_MDSE"]["MDSE_T"]) # MDSE_T at optimum print(opt["y"]["anc"]["opt_MDSE"]["MDSE_S"]) # MDSE_S at optimum # ################################################# # Table 2, Column 3 # # ################################################# # Random saturation design (two partially treated clusters), ICC = 0.0 # Optimal cluster saturations and allocation n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 1.0 # variance of cluster error (tau^2) sigma = 0.0 # variance of individual error (sigma^2) # BRUTE FORCE SEARCH ##x = np.arange(0.01,.5+.005,.01) ##MDSE_T = np.zeros((len(x),len(x))) ##MDSE_S = copy.deepcopy(MDSE_T) ##for i in range(len(x)): ## pi = [x[i],1-x[i]] ## for j in range(len(x)): ## f = [x[j],1-x[j]] ## slope = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) ## X = [slope[k] for k in ["MDSE_T","MDSE_S"]] ## MDSE_T[i,j] = copy.deepcopy(X[0]) ## MDSE_S[i,j] = copy.deepcopy(X[1]) ## ## ##theta = 1/2 # fraction of weight put onto MDSE_T ##obj = theta*MDSE_T + (1-theta)*MDSE_S; # objective function to maximize ##ind = index_min(obj); # index of psi that minimizes objective ##print(x[ind[0]]) # optimal saturation for lower saturation ##print(1 - x[ind[0]]) # optimal saturation for higher saturation ##print(x[ind[1]]) # optimal fraction of clusters in lower saturation ##print(1 - x[ind[1]]) # optimal fraction of clusters in higher saturation ##print(MDSE_T[ind[0], ind[1]]) # MDSE_T at optimum ##print(MDSE_S[ind[0], ind[1]]) # MDSE_S at optimum # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0] < .5 and 0 < x[1] and x[1] < .5 def objective(x): # x contains the first saturation level, and the number of people in the least saturated group. pi = [x[0],1-x[0]] f = [x[1],1-x[1]] value = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,1) objective_function = theta*value["MDSE_T"] + theta*value["MDSE_S"] return {"obj":objective_function, "anc":{"opt_MDSE":value}} opt = Compass_Search(objective, admissable, x0=np.array([.25,.25])) print(opt["x"][0]) # optimal saturation for lower saturation print(1-opt["x"][0]) # optimal saturation for higher saturation print(opt["x"][1]) # optimal fraction of clusters in lower saturation print(1-opt["x"][1]) # optimal fraction of clusters in higher saturation print(opt["y"]["anc"]["opt_MDSE"]["MDSE_T"]) # MDSE_T at optimum print(opt["y"]["anc"]["opt_MDSE"]["MDSE_S"]) # MDSE_S at optimum # ################################################# # Table 2, Column 4 # # ################################################# # Joint maximization across MDE and MDSE. # NOTE: pi has been maximized using the a joint maximization process, # performed in Complete_Maximization.m. n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) # BRUTE FORCE SEARCH ##pi = [0,0.21,0.88] # cluster saturations. ##x = np.arange(0.01,.99,.01) # values for maximization. ##MDE_T= np.ones((len(x), len(x)))*np.inf ##MDE_S = copy.deepcopy(MDE_T) ##MDSE_T = copy.deepcopy(MDE_T) ##MDSE_S = copy.deepcopy(MDE_T) ## ##for i in range(len(x)): ## for j in range(len(x)-i): ## f = copy.deepcopy([x[i],x[j],1-x[i]-x[j]]) ## pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## MDE_T[i,j] = copy.deepcopy(pooled["MDE_T"]) ## MDE_S[i,j] = copy.deepcopy(pooled["MDE_S"]) ## MDSE_T[i,j] = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,1,2)["MDSE_T"] ## MDSE_S[i,j] = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,2)["MDSE_S"] ## ##obj = MDE_T+MDE_S+MDSE_T+MDSE_S; # objective function to maximize ##ind = index_min(obj); # indeces of x that minimizes objective function ##print([x[ind[0]], x[ind[1]], 1-x[ind[0]]-x[ind[1]]]) # optimal pi ##print([MDE_T[ind[0], ind[1]],MDE_S[ind[0], ind[1]],MDSE_T[ind[0], ind[1]],MDSE_S[ind[0], ind[1]]]) # optimal MDEs/MDSEs. # COMPASS SEARCH def admissable(x): return 0 < x[0] and x[0] < x[1] and x[1] < 1 and 0 < x[2] and 0 < x[3] and x[2]+x[3] < 1 def objective(x):# x contains the second two saturation levels, and the number of people in the first two control groups. pi = [0, x[0], x[1]] f = [x[2], x[3], 1-x[2]-x[3]] pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) slope1 = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,1,2)["MDSE_T"] slope2 = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,2)["MDSE_S"] objective_function = slope1 + slope2 + pooled["MDE_T"] + pooled["MDE_S"] return {"obj":objective_function, "anc":{"pooled":pooled, "slope1":slope1, "slope2":slope2}} opt = Compass_Search(objective, admissable, x0=np.array([.25,.5, .5, .25])) print(0, *opt["x"][:2]) # optimal pi print(opt["x"][2],opt["x"][3],1-sum(opt["x"][2:])) # optimal f print(opt["y"]["anc"]["pooled"]["MDE_T"], opt["y"]["anc"]["pooled"]["MDE_S"], opt["y"]["anc"]["slope1"], opt["y"]["anc"]["slope2"]) # ################################################# # Table 2, Column 5 # # ################################################# # MDEs resulting from study design of Banerjee et al. n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) pi = [0, 0.25, 0.5, 0.75, 1] # cluster saturations f = [0.2, 0.2, 0.2, 0.2, 0.2] # fraction of clusters in each saturation MDE_T = 0; MDE_S = 0; MDE_Tonly = 0; MDSE_T = 0; MDSE_S = 0 pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) slope_T = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,1,3) slope_S = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,0,2) print(pooled["MDE_T"], pooled["MDE_S"], slope_T["MDSE_T"], slope_S["MDSE_S"]) # ################################################# # Table 2, Column 6 # # ################################################# # MDEs resulting from study design of Sinclair et al. n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) pi = [0, 0.1, 0.5, 1] # cluster saturations f = [0.25, 0.25, 0.25, 0.25] # fraction of clusters in each saturation MDE_T = 0; MDE_S = 0; MDE_Tonly = 0; MDSE_T = 0; MDSE_S = 0 pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) slope_T = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,1,3) slope_S = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,0,2) print(pooled["MDE_T"], pooled["MDE_S"], slope_T["MDSE_T"], slope_S["MDSE_S"]) # ################################################# # Table 2, Column 7-8 # # ################################################# # MDEs resulting from study design of Baird et al. n = 10 # cluster size C = 100 # number of clusters alpha = 0.05 # significance gamma = 0.8 # power tau = 0.1 # variance of cluster error (tau^2) sigma = 0.9 # variance of individual error (sigma^2) pi = [0, 0.33, 0.67, 1] # cluster saturations f = [0.55, 0.15, 0.15, 0.15] # fraction of clusters in each saturation MDE_T = 0; MDE_S = 0; MDE_Tonly = 0; MDSE_T = 0; MDSE_S = 0 pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) slope_T = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,1,3) slope_S = power_slope (n,C,alpha,gamma,tau,sigma,pi,f,0,2) print(pooled["MDE_T"], pooled["MDE_S"], slope_T["MDSE_T"], slope_S["MDSE_S"]) #### Table 2, Column 8 # Design with optimal MDEs subject to: # 1.) MDE_T held no larger than that of Baird et al. (MDE_t <= 0.2679) # 2.) Saturation bins held the same as Baird et al. # 3.) The share of cluster in two middle saturations held equal. # BRUTE FORCE SEARCH ##x = np.arange(0.01, .99+.005, .01) # values for optimization ##MDE_T= np.ones((len(x), len(x)))*np.inf ##MDE_S = copy.deepcopy(MDE_T) ##MDSE_T = copy.deepcopy(MDE_T) ##MDSE_S = copy.deepcopy(MDE_T) ## ##for i in range(len(x)): ## for j in range(len(x)-i): ## f = [x[i],(1-x[i]-x[j])/2,(1-x[i]-x[j])/2, x[j]] ## X = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) ## MDE_T[i,j] = X["MDE_T"]; MDE_S[i,j] = X["MDE_S"] ## MDSE_T[i,j] = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,1,3)["MDSE_T"] ## MDSE_S[i,j] = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,2)["MDSE_S"] ## ##MDE_S_aug = copy.deepcopy(MDE_S) ##for i in range(len(x)): ## for j in range(len(x)): ## if MDE_T[i,j] > 0.2679: ## MDE_S_aug[i,j] = np.inf ##obj = MDE_S_aug # objective function to maximize ##ind = index_min(obj) # index of x that minimizes MDE_S ##print(x[ind[0]]) # optimal fraction of clusters in pure control ##print((1-x[ind[0]]-x[ind[1]])/2) # optimal fraction of clusters in two middle saturations ##print(x[ind[1]]) # optimal fraction of clusters in pure treatment ##print(MDE_T[ind[0], ind[1]]) # MDE_T at optimal control size ##print(MDE_S[ind[0], ind[1]]) # MDE_S at optimal control size ##print(MDSE_T[ind[0], ind[1]]) # MDSE_T at optimal control size ##print(MDSE_S[ind[0], ind[1]]) # MDSE_S at optimal control size # COMPASS SEARCH def minimum_admissable(x): return 0 < x and x < 1 def minimum_objective(theta): def admissable(x): return 0 < x[0] and 0 < x[1] and x[0] + x[1] < 1 def objective(x):# x contains the first and last saturation level. f = [x[0], (1-x[0]-x[1])/2,(1-x[0]-x[1])/2, x[1]] pooled = power_pooled(n,C,alpha,gamma,tau,sigma,pi,f) slope1 = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,1,3)["MDSE_T"] # MDSE_T at optimal control size slope2 = power_slope(n,C,alpha,gamma,tau,sigma,pi,f,0,2)["MDSE_S"] # MDSE_S at optimal control size objective_function = theta*pooled["MDE_T"] + (1-theta)*pooled["MDE_S"] return {"obj":objective_function, "anc":{"x":x, "pooled":pooled, "slope1":slope1, "slope2":slope2}} opt = Compass_Search(objective, admissable, x0 = np.array([0.25, 0.25])) return {"obj":abs(opt["y"]["anc"]["pooled"]["MDE_T"]-0.2679), "anc":{"opt":opt}} theta = Compass_Search(minimum_objective, minimum_admissable, x0=np.array([1])) print(theta["y"]["anc"]["opt"]["y"]["anc"]["x"][0]) # optimal fraction of clusters in pure control x=theta["y"]["anc"]["opt"]["y"]["anc"]["x"] print((1-x[0]-x[1])/2) # optimal fraction of clusters in two middle saturations print(theta["y"]["anc"]["opt"]["y"]["anc"]["x"][1]) # optimal fraction of clusters in pure treatment print(theta["y"]["anc"]["opt"]["y"]["anc"]["pooled"]["MDE_T"]) # MDE_T at optimal control size print(theta["y"]["anc"]["opt"]["y"]["anc"]["pooled"]["MDE_S"]) # MDE_S at optimal control size print(theta["y"]["anc"]["opt"]["y"]["anc"]["slope1"]) # MDSE_T at optimal control size print(theta["y"]["anc"]["opt"]["y"]["anc"]["slope2"]) # MDSE_S at optimal control size